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ABSTRACT 


We report the discovery of a sextuply-eclipsing sextuple star system from T'ESS data, TIC 168789840, 
also known as TYC 7037-89-1, the first known sextuple system consisting of three eclipsing binaries. 
The target was observed in Sectors 4 and 5 during Cycle 1, with lightcurves extracted from TESS Full 
Frame Image data. It was also previously observed by the WASP survey and ASAS-SN. The system 
consists of three gravitationally-bound eclipsing binaries in a hierarchical structure of an inner quadru- 
ple system with an outer binary subsystem. Follow-up observations from several different observatories 
were conducted as a means of determining additional parameters. The system was resolved by speckle 
interferometry with a 0742 separation between the inner quadruple and outer binary, inferring an esti- 
mated outer period of ~2 kyr. It was determined that the fainter of the two resolved components is an 
8.217 day eclipsing binary, which orbits the inner quadruple that contains two eclipsing binaries with 
periods of 1.570 days and 1.306 days. MCMC analysis of the stellar parameters has shown that the 
three binaries of TIC 168789840 are “triplets”, as each binary is quite similar to the others in terms 
of mass, radius, and Tę. As a consequence of its rare composition, structure, and orientation, this 
object can provide important new insight into the formation, dynamics, and evolution of multiple star 
systems. Future observations could reveal if the intermediate and outer orbital planes are all aligned 
with the planes of the three inner eclipsing binaries. 


Keywords: Eclipsing Binary Stars — Transit photometry — Astronomy data analysis — Multiple star 


systems — Machine learning — High-performance computing 


1. INTRODUCTION 


The Transiting Exoplanet Survey Satellite (TESS) 
mission (Ricker et al. 2015) has dramatically improved 
our ability to discover multiple star systems. Though it 
is more prone to systematics than the Kepler telescope 
and has a poorer angular resolution (21" per pixel for 
TESS vs. 3.98" per pixel for Kepler), the breadth of ob- 
servation of T'ESS, encompassing nearly the entire sky, 
has allowed for the identification of many candidate mul- 
tiple star systems through the analysis of eclipses in the 
lightcurves. In fact, a collaboration between the NASA 
Goddard Space Flight Center (GSFC) Astrophysics Sci- 
ence Division and the MIT Kavli Institute, in conjunc- 
tion with expert visual surveyors, has found well over 
100 triple and quadruple star system candidates. This 
number will continue to increase as T'ESS proceeds with 
the extended mission at faster observation cadence (10 
minutes for cycles 3 and 4 vs. 30 minutes for cycles 1 
and 2), enabling researchers to capture shorter-duration 
eclipse events. We also note that lightcurves from cycles 
1 and 2 have yet to be fully exploited. 


The large majority of our discovered candidate triple 
and quadruple star systems are quadruples, followed by 
triples. Though quadruple systems are much more rare 
than triple systems, the large outer orbit of the third 
star in a hierarchical triple, necessary for stability, sub- 
stantially reduces the probability that the eclipse or oc- 
cultation of the third star will be visually noticed in a 
TESS lightcurve. Beyond quadruple stars, the proba- 
bility of systems with more components being identified 
via photometry alone is remote, as the formation of sex- 
tuple systems is likely quite rare. This low probabil- 
ity is compounded by the requirement that each binary 
must be oriented in such a manner that they are all 
eclipsing. Though simulations of stellar system forma- 
tion have found that a sextuple system consisting of two 
inner triples is nearly ten times more likely to form than 
a system of three close binaries (van den Berk et al. 
2007), the visual detection of all the eclipses in a sextu- 
ple consisting of two triples is far less likely, again, due 
to the large outer orbit of the third star in each triple. 

In this work we present a sextuple system which ex- 
hibits all six eclipses (three primary and three sec- 
ondary) discovered with T'ESS. We show that TIC 


168789840 consists of three close binaries. The inner 
quadruple system with a period of ~3.7 yrs is comprised 
of two eclipsing binaries (which we provide the names 
“A” and “C”), at periods of 1.570 and 1.306 days, re- 
spectively; the inner quadruple is orbited by another 
eclipsing binary (which we call ^B"), with a period of 
8.217 days, at a period of ~2 kyr. The structure of 
the system, shown in Figure 1, will be the nomenclature 
that will be used for the rest of this paper. Prior to the 
discovery of TIC 168789840, there were 17 known sex- 
tuple star systems according to the June 2020 update 
of the Multiple Star Catalog (Tokovinin 2018a). TIC 
168789840 is the first that is sextuply eclipsing, with 
the caveat that Jayaraman, Rappaport, Borkovits, Za- 
sche et al. are currently analyzing another such system 
that will be published in the near future. 

There are several known sextuple systems with a sim- 
ilar structure to that of TIC 168789840. One of these 
systems, ADS 9731, is a resolved visual quadruple sys- 
tem known for more than a century; two of its compo- 
nents were further determined to be close spectroscopic 
binaries by Tokovinin et al. (1998). 88 Tauri, suspected 
to be a spectroscopic quintuple by Burkhart & Coupry 
(1988), was later determined to be a sextuple system of 
three binaries (Tokovinin 1997), with the two binaries 
comprising the inner quadruple having an 18 yr period 
(Lane et al. 2007). Interestingly, of the known sextuple 
systems, TIC 168789840 is most similar to the famous 
Castor system, which also contains three close binaries. 

Castor, among the brightest star systems in the sky, 
was originally identified as a visual binary system in 
1719 by Bradley and Pound. Belopolsky (1897) found 
that one of its components was a spectroscopic binary, 
and Curtis (1905) discovered that another component 
was also a binary. Adams & Joy (1920) found that there 
was a third component which was also a binary, complet- 
ing the discovery of the Castor system as the first known 
sextuple star system. The mass and radius ratios of the 
binaries of TIC 168789840, in addition to the close or- 
bits of the binaries, are found in this work to be quite 
similar to those determined by the extensive analysis of 
Castor. 

The rest of the paper is organized as follows. Sec- 
tion 2 outlines the initial detection and analysis of the 
TESS data; the disentaglement of the individual EB 
lightcurves is presented in Section 3. In Sections 4 and 5 
we present the analysis of archival data and our follow- 
up observations, respectively. The comprehensive analy- 
sis of the parameters of the system and the correspond- 
ing discussion of the results is presented in Section 6. 
Finally, we draw our conclusions in Section 7. 


2. DETECTION 


Using the 129,000-core Discover supercomputer at 
the NASA Center for Climate Simulation (NCCS) at 
NASA GSFC, we are building Full-Frame-Image (FFI) 
lightcurves for all stars observed by TESS up to 15th 
magnitude. All original and calibrated FFIs are pro- 
duced by the TESS Science Processing Operations 
Center (SPOC, Jenkins et al. 2016). Target lists 
were created through a parallelized implementation of 
tess-point (Burke et al. 2020) on the TESS Input Cat- 
alog (TIC) provided by the Mikulski Archive for Space 
Telescopes (MAST). The lightcurves for each sector were 
constructed in 1-4 days of wall clock time (for a total of 
over 100 CPU-years to date), depending on the den- 
sity of targets in the sector, through a parallelized im- 
plementation of the eleanor Python module (Feinstein 
et al. 2019). From these lightcurves, we are performing 
a search for multiple stellar systems using targets from 
the GSFC TESS Eclipsing Binary (EB) Catalog (Kruse 
et al. 2021, in prep). 

This catalog of eclipsing binaries was generated by 
a neural network classifier. This neural network was 
trained on the NCCS Advanced Data Analytics PlaT- 
form (ADAPT) GPU cluster to classify a lightcurve 
(as either an EB or not an EB) based only on 
the feature of the eclipse, neglecting any periodic- 
ity or time-dependency. The neural network is a 
one-dimensional adaptation of the ResNet (He et al. 
2015) structure to accommodate the data shape of a 
lightcurve, built in Python using keras (Chollet et al. 
2015)/tensorflow(Abadi et al. 2015). A strength of 
this approach is that it allows for the identification of 
single-eclipse EBs. As such, a lightcurve with an eclipse 
recognizable by the neural network, no matter the num- 
ber of eclipses occurring in a single lightcurve, will be 
properly classified as an EB by the neural network. Fig- 
ure 2 shows the activation of the neural network on the 
feature of the eclipse in a segment of the TIC 168789840 
lightcurve, which demonstrates that each eclipse does 
not need to be individually identified by the neural net- 
work in order for the lightcurve to be classified as an EB. 
The lack of a periodicity or similarity constraint allows 
for a lightcurve with multiple irregular eclipses, such as 
TIC 168789840, to be classified as an EB. 

Lightcurves with multiple sets of eclipses are manu- 
ally flagged as meriting further investigation. While the 
overwhelming majority of these lightcurves are deter- 
mined to be false positives caused by close proximity of 
two or more EBs blending into a single lightcurve, there 
remains a fraction which cannot be explained by such 
contamination. This is determined through photocen- 
ter analysis, the output of which for TIC 168789840 is 
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Figure 1. Structure of TIC 168789840, a sextuple system of three eclipsing binaries arranged as inner quadruple AC and outer 
binary B. In this work we will discuss how we arrived at this configuration. 
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Figure 2. Saliency map - indicating features of importance to the classification — of a segment of the TIC 168789840 eleanor 
raw flux lightcurve. The neural network activates on the feature of the eclipse. Some eclipses activate more strongly than 
others, which is a function of the response of the neural network to the lightcurve as well as idiosyncrasies of the training data. 
The map was created using keras-vis (Kotikalapudi & contributors 2017) on the penultimate layer of the EB classifier neural 


network. 


shown in Figure 3. The analysis follows the difference 
imaging procedure of Bryson et al. (2013) as adapted 
into the DAVE vetting pipeline (Kostov et al. 2019). 
Briefly, we first perform a Box Least Squares (BLS) 
analysis (Kovacs et al. 2002) of the TESS lightcurve 
to measure the ephemerides for all sets of eclipses. As 
demonstrated in Figure 4, the TESS lightcurve of TIC 
168789840 shows three distinct periods with primary 
and secondary eclipses. Next, for each eclipse of each 
set we create the out-of-eclipse images and difference im- 


ages, measure the corresponding center-of-light (photo- 
center) by calculating the respective x- and y-moments 
and, finally, compare the average out-of-eclipse photo- 
center to the average difference image photocenter for 
each set. A significant shift between these two photo- 
centers would indicate a potential false positive due to a 
nearby field star. For TIC 168789840, there are no sig- 
nificant differences between the measured out-of-eclipse 


Figure 3.  Photocenter analysis for the three primary 
eclipses of TIC 168789840 for Sector 4. The panels show 
the mean difference image for pair A (top), B (middle) and 
C (bottom). The large open circle represents the average dif- 
ference image photocenter and the large star symbol repre- 
sents the catalog position of the target. The small x symbols 
represent nearby stars from the TIC within 10 magnitude 
difference. The difference image photocenters for all three 
sets of eclipses are on-target. 


and difference-image photocenters for all sets of eclipses, 
indicating that the target is their source. 


3. DISENTANGLEMENT OF THE LIGHTCURVES 
3.1. Fourier Method 


5 


After identifying the sources of all the eclipses to 
be on target (i.e. belonging to TIC 168789840), we 
needed to disentangle the combined photometry to cre- 
ate lightcurves for each of the three eclipsing binaries. 
Here we introduce one of our two methods for disen- 
tangling the photometric lightcurves of the three eclips- 
ing binaries, superposed in the TESS data. This ap- 
proach, which amounts to a Fourier decomposition, re- 
quires prior knowledge of the orbital periods. In this 
particular case, the periods were determined from Lomb- 
Scargle and BLS transforms of the TESS data as well as 
archival ASAS-SN data (see Section 4.1). 

To represent the three binaries in the lightcurve, we 
fit a harmonic series of the following form to the entire 
27-day TESS data train: 


3 / 50 
F(t) = ` p o) sin(u,t) + Bt? «ou +7, 


n=1 

(1) 
where wp is the nth orbital frequency in the series repre- 
senting the mth binary, and is given by 27n/P,,, where 
Pm is the orbital period of the mth binary. In all there 
are 3 x 50x 2+ 1 = 301 linear coefficients to be fit, 
i.e., all the o, 8,, and y (the latter being the constant 
background level, which is ~ 1 if the lightcurve is nor- 
malized). We note that the values of wn are unrelated to 
the usual orthogonal frequencies used in an FFT which 
are given by integer multiples of 27/T, where T is the 
duration of the observation interval. 

'These coefficients can all be fitted simultaneously with 
the inversion of a single 301 x 301 y? matrix, which takes 
much less than a minute on a standard laptop. While we 
used 50 harmonics in this case, we have found that 30 
harmonic terms are sufficient to effectively reconstruct 
most binary lightcurves, except those with very deep 
and/or sharp eclipses. 

'The next and final step in the procedure is to recon- 
struct the lightcurve for the mth binary via the following 
sum: 


50 
Fm(tj) = b» al™ sin(wntj) + BL cos(wst;) +y (2) 
n=1 


where 7 is the 7th data point. 

The results of the Fourier disentanglement for the 
TESS lightcurve of TIC 168789840 are shown in Fig- 
ure 5. The three panels show the reconstructed 
lightcurves for the A, B, and C binaries with periods 
of 1.570 days, 8.217 days, and 1.306 days, respectively. 
These are perfectly conventional EBs with two eclipses 
per period, and, at first glance, the lightcurves seem to 
indicate circular orbits. 
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Figure 4. Upper panel: The TESS lightcurve of TIC 168789840 for sectors 4 and 5. Several eclipses are blended, most notably 
around time 1461. The rest of the panels show the lightcurve phase-folded on the three distinct periods of the three EBs as 


measured by BLS. 


Finally, we discuss an important caveat to this 
Fourier-based method for disentangling multiple super- 
posed lightcurves. This technique works best if none of 
the harmonics of one eclipsing binary overlaps, within a 
resolution element (2x/T), of any of the harmonics from 
the other binaries. If there is significant overlap among 
any of the lower harmonics (e.g., for n < 5) then this 
technique may have problems with the degenerate fre- 
quencies. If the overlap is among the higher harmonics 


(e.g., for n Z 15), then this effect is probably negligible. 
In the case of TIC 168789840 the 4th harmonic of the 
A binary overlaps the 21st harmonic of B, while the 5th 
harmonic of the C binary overlaps the 6th harmonic of 
binary A. We have checked what problems this might 
cause, by removing each of the binaries separately, and 
we find very similar results to fitting for them all simul- 
taneously. 


3.2. Iterative Method 
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Figure 5. Reconstructed T'ESS lightcurves for the A, B, and C binaries. These are presented on the same y-axis in order to 


visualize the relative contributions from each binary. 


We have also used another independent method for 
disentangling the lightcurves. The results of the two 
procedures can be used to check each other. 

In the iterative approach, we follow the schematic 
steps outlined in Table 1. We start with the original 
lightcurve time series denoted as TABC where the “T” 
stands for time series, and the ‘A’, ‘B’, and ‘C’ signify 
that all three binaries are contained in T. We then start 
by producing a phase-folded, binned and averaged or- 
bital lightcurve (hereafter, denoted as ‘F’) for one of 
the binaries (e.g., A), by first removing from the time 
series the intervals when eclipses from the other two bi- 
naries (e.g., B and C) occur. The residual time series 
is then phase folded, binned and averaged to produce 
FA, ;. The subscripts ‘a’ and ‘T’ signify that we are 
starting the cleaning process by working along track ‘a’ 
(see Table 1 for the track definitions) and this will be a 
first-level product (‘T’). 

The next step is to subtract fold FA, ; from TABC, us- 
ing a three-point local Lagrangian-interpolation to cal- 
culate the flux to be subtracted at each observed pho- 
tometric phase of the binary. The result is a time series 
comprised of the blended lightcurves of only two of the 
three binaries, and we denote this product as, e.g., TBC. 
This completes the first-level products. In all there are 
three preliminary folds, one for each binary, and three 
time series!, each containing two of the binaries. See 
the second row of Table 1. 

To produce the second-level (‘II’) products, we take 
each of the time series from level I, comprised of two bi- 
naries, e.g., TBC, remove the eclipses of either one of the 
binaries, and produce a phase-folded, binned, and aver- 
aged orbital lightcurve of the other binary, e.g., FB, rr. 
Then, as in level I, we subtract off the folded lightcurve 
of that binary to produce a time series containing only a 


l Note, for practical reasons, we added a constant flux to these 
time series in such a way that the flux of the very first data point 
retained the same value as in the original time series. In this 
manner, we replaced the varying light of the extracted binary 
with a constant extra light. 


single binary. Schematically, TBC - FB, ;; = TC, and 
TBC - FC,,;; = TB. The net result of the level II prod- 
ucts are two semi-independent folded orbital lightcurves 
for the A, B, and C binaries (six folds in all), and two 
semi-independent time series for binaries A, B, and C 
(six in all). See the third row of Table 1 for the full set 
of second-level products. 

'The final step is to take all six time series and fold 
them about the orbital period of the single binary re- 
maining in each one. This yields two semi-independent 
pairs of phase folded orbital lightcurves, e.g., FC, rr; 
and FC, zzz for each binary. Refer to the fourth row of 
Table 1 for the set of final folded lightcurves. 

We applied the complete iterative disentangling 
method to two different initial time-series. The first was 
for the original time-series obtained from the TESS data 
with the use of the Fitsh pipeline of Andras Pál (Pal 
2012). Second, in order to reduce the non-physical scat- 
ter of the extracted lightcurves, we removed a 6-day-long 
section of the lightcurve between BJD 2458418.4 and 
2458424.7 due to its large slope and, furthermore, we 
carried out a minor detrending operation with the soft- 
ware package of Wotan (Hippke et al. 2019) to remove 
some additional, slight flux-level variations on a time- 
scale of 10-15 days. In this manner, the noise-level of 
the disentangled lightcurves was reduced significantly, 
without any changes in the structure. "Therefore, for 
our analysis, we used the data series obtained from this 
slightly detrended second time-series. 


4. ARCHIVAL DATA 


A search for archival data on TIC 168789840 reveals 
that there are a couple of rich sources of historical pho- 
tometry. Figure 6 highlights the baseline covered by the 
available archival observations of the target from ASAS- 
SN, WASP and TESS; the corresponding ephemerides 
for the three EBs are listed in Table 2. 


4.1. All-Sky Automated Survey for Supernovae 
(ASAS-SN) 


Table 1. Logical Tree for Iterative Disentanglement 
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Figure 6. Archival data from ASAS-SN (blue and green symbols), WASP 200-mm lenses (magenta symbols), WASP 85-mm 
lenses (turquoise symbols), and T'ESS (red symbols) highlighting the baseline covered by the photometry. 


TIC 168789840 was observed by the ASAS-SN 
(Shappee et al. 2014, Kochanek et al. 2017) Sky Patrol 
from BJD 2456000 to 2459100, with excellent coverage 
over the last 7 years. In all, there are 4746 archival pho- 
tometric data points available. After renormalizing the 
green-band to the visual-band observations, we carried 
out a BLS (Kovács et al. 2002) transform of the data 
to see which of the three binary EBs we could recover. 
The top two highest peaks in the BLS transform were 
of the 1.570 day and 8.217 day periods (from binaries A 
and B, respectively). We then used the Fourier cleaning 
tool described in Section 3.1 to remove these two peri- 
ods from the data. The BLS transform of the cleaned 
ASAS-SN lightcurve then reveals the 1.306 day primary 
eclipses of the C binary. In the top panels of Figure 7 we 
show the ASAS-SN data folded about the three periods 
determined from these data. 


4.2. Wide-Angle Search for Planets (WASP) 


The field of TIC 168789840 was observed by the 
WASP-South transit search between 2006 and 2014. 
WASP-South was an array of 8 cameras located in 
Sutherland, South Africa (Pollacco et al. 2006). Be- 
tween 2006 and 2011 the cameras used 200-mm, f/1.8 
lenses, observing with a 400—700-nm filter and using a 
48" photometric extraction aperture. Between 2012 and 
2014 the cameras had 85-mm, f/1.2 lenses with an SDSS- 
r filter and a 112" extraction aperture. TIC 168789840 
is the brightest star in both-size apertures (the next 


brightest is 2.5 magnitudes fainter). Observations had 
a typical 12-min cadence, and where obtained on clear 
nights spanning 150 days in each of 2006, 2007, 2011, 
2012, 2013 and 2014. A total of 126000 photometric 
data points were recorded. However, we found that the 
S/N was better using only the 18,000 data points taken 
with the 200-mm lens. 

We analyzed the WASP data in the same manner that 
we did for the ASAS-SN data. Again, the eclipses from 
the A and B binaries were the easiest to find. We then 
cleaned the data of these two periods, and easily de- 
tected the eclipses of the C binary. The bottom panels 
of Figure 7 show the WASP data folded in the same 
manner as the ASAS-SN data. 


5. FOLLOW-UP OBSERVATIONS 


Upon identification of the system, we had overwhelm- 
ing support from follow-up observers providing nearly 
fifty separate measurements from seven different obser- 
vatories. These range from photometric measurements 
to radial velocity and speckle imaging, each helping us 
to further unravel the nature of the system. 


5.1. Photometric measurements 
5.1.1. TESS Followup Observing Program 
Photometric follow-up observations were performed 
through Subgroup 1 of the TESS Follow Up Observing 


Program (TFOP) as described in more detail below. We 
used the TESS Transit Finder, which is a customized 
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Figure 7. Folds of the ASAS-SN (top) and WASP (bottom) archival data for the A, B, and C binaries about their respective 
orbital periods. For the A and B folds, the raw data were used. Before constructing the fold about the C period we removed 
the orbital profiles of the A and B binaries by Fourier filtering (see section 3.1). For each binary we have written the fold period 
and the epoch of the primary eclipse on the plot. The fold has been shifted by half an orbital period for aesthetic reasons. 


Table 2. Derived ephemerides for the three EBs in TIC 168789840 from ASAS-SN, TESS and WASP. 


Binary A B C 
TESS 
Period [days] 1.570101 8.217173 1.305934 
TO [BJD - 2450000] 8412.3855 8411.9008 8413.6822 
ASAS — SN 
Period [days] 1.569984 8.216958 1.305904 
TO [BJD - 2450000] 5950.642 5946.786 5946.843 
WASP 
Period [days] 1.570044 8.217670 1.305878 
TO [BJD - 2450000] 3900.150 3900.699 3900.564 
Radial Velocities fixed fixed fixed 
TO [BJD - 2450000] 9151.868 9151.446 9151.193 
Global Fitted Periods^ | 1.570013(9) | 8.217111(30) | 1.305883(6) 


Notes. (a) The long-term average period is determined from a linear fit to the four independently determined times of eclipse. 
In the case of binary B this assumes no change in its center-of-mass velocity over the past 15 years. In the case of binaries A 
and C, which we later show to be in a ~3.7 year quadruple orbit, with speeds of ~7 km s^ !, this could lead to effects as large 

as 23 parts per million in the reported period. But, much of the latter is averaged over in the WASP and ASAS-SN 
measurements which span the ~3.7 year orbit. 


version of the Tapir software package (Jensen 2013), to 
schedule our transit observations. These observations, 
shown in Figure 8, confirm that the target is the source 
of the different sets of eclipses detected in TESS data 
and rule out contamination from nearby sources. Several 
of the observations shown in Figure 8, while targeted at 
one particular eclipse of a given binary, simultaneously 
observed eclipses from either of the other two binaries. 


'TIC 168789840 was observed on nine nights with the 


Evans telescope at El Sauce in Coquimbo Province, 
Chile. This system consists of a 0.36-m CDK telescope 
with a SBIG STT1603-3 CCD, which has an image scale 
of 1747 pixel! and 18/8 x 12/5 field of view; all obser- 
vations used the Rc filter. The observations covered 
the following UTC dates and eclipses: 2020-10-07 (C 
primary); 2020-10-18 (A primary); 2020-10-19 (A sec- 
ondary, C secondary); 2020-10-21 (C primary); 2020-10- 
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Figure 8. TFOP-led photometric observations of TIC 
168789840 confirming that the target is the source of the 
eclipses detected in T'ESS data. Primary eclipses of A, B, 
and C are shown in the three panels. The vertical red bands 
represent the corresponding ingress and egress times. The in- 
dividual measurements are vertically offset for clarity. Some 
observations cover two eclipses (see text). 


22 (A secondary); 2020-10-23 (A primary, C secondary); 
2020-11-02 (A secondary); 2020-11-06 (A primary, B pri- 
mary); and 2020-11-10 (A secondary, B secondary). The 
photometric data were extracted using the AstroImageJ 
(AIJ) software package (Collins et al. 2017). 

The Perth Exoplanet Survey Telescope (PEST) is a 
0.3 m telescope in Perth, Australia, with an image scale 
of 172 and a 31’ x 21' field of view. PEST observed 
in the Ro filter on UTC 2020-10-20, covering the B 
primary and A secondary eclipses. A custom pipeline 
based on C-Munipack (Motl 2011) was used to calibrate 
the images and extract the differential photometry. 

Two observations made use of the Las Cumbres Ob- 
servatory Global Telescope (LCOGT) network (Brown 
et al. 2013). Primary eclipses of both the A and C bi- 
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Figure 9. Phase folded lightcurve of ground-based FRAM 
data showing binary A (filter R was used) plotted against 
the PHOEBE fit. 


naries were observed on UTC 2020-10-19 using a 0.4- 
m telescope in Sutherland, South Africa. The LCOGT 
0.4-m telescopes are equipped with 2048 x 3072 SBIG 
STX6303 cameras having an image scale of 0/57 pixel~! 
resulting in a 19’ x 29’ field of view. On 2020-11-08, one 
of the LCOGT 1.0-m telescopes at Siding Spring Ob- 
servatory observed this system, covering the C primary 
and A secondary eclipses. The 4096 x 4096 LCOGT 
SINISTRO cameras have an image scale of 07389 per 
pixel, resulting in a 26’ x 26’ field of view. The LCOGT 
images were calibrated by the standard LCOGT BAN- 
ZAI pipeline (McCully et al. 2018) and the photometric 
data were extracted using AIJ. 


5.1.2. FRAM 


Some follow-up photometric data from ground based 
observatories were also obtained with a small 30-cm 
telescope FRAM. It is the Orion ODK 300/2040mm, 
equipped with the CCD camera MII G4-16000. All ob- 
servations carried out in standard R filter. The FRAM 
telescope itself (Janecek et al. 2019) is located at the 
peak of Los Leones, near the town of Malargüe, at the 
Pierre Auger Observatory, Argentina. The phase fold of 
ten separate observations is shown in Figure 9. 


5.2. Spectroscopy 
5.2.1. CHIRON 


Eight high-resolution optical spectra of TIC 
168789840 were taken with the CHIRON fiber echelle 
spectrometer (Tokovinin et al. 2013) at the CTIO 1.5 m 
telescope operated by the Small & Moderate Aperture 
Research Telescope System (SMARTS) consortium be- 
tween 2020 November 6 and 21. The spectral resolution 
is 80,000, exposure time 15 min., and the typical S/N 
is ~20 per pixel (pixel width ~1.2 km s^). The wave- 
length calibration is determined from the ThAr spectra 
taken immediately after the stellar spectra. 
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Figure 10. Top panel: Cross-correlation functions (CCFs) 
of the 8 CHIRON spectra with a binary mask. Bottom panel: 
Same for the 4 TRES spectra. The Julian dates of observa- 
tions are indicated. The RVs of the components correspond- 
ing to their orbits are marked by the colored ticks above 
each CCF. The dip at zero RV in the first T'RES spectrum 
is produced by Moon contamination. 


To find the radial velocities (RVs), the spectra were 
cross-correlated with a binary mask based on the so- 
lar spectrum. Details of this procedure are provided 
in (Tokovinin 2016). Only wavelengths from 480nm 
to 650nm are used. The cross-correlation functions 
(CCFs) show a narrow dip with an amplitude of 0.07- 
0.08 and an rms width of 7 km s^! that corresponds to 
the projected rotation speed of 10.2 km s-!, see Fig- 
ure 10. Moreover, there are broad features resulting 
from other components with fast rotation. Analysis of 
CHIRON spectra shows convincingly that the narrow 
dip belongs to the primary component of the 8-day pair, 
B1. A circular spectroscopic orbit fits both CHIRON 
and TRES RVs. The RVs of the rapid rotators A1 and 
C1 are derived by modeling the spectrum (Section 5.3). 
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5.2.2. Tillinghast Reflector Echelle Spectrograph (TRES) 


Four additional spectroscopic observations were made 
with the Tillinghast Reflector Echelle Spectrograph 
(TRES; Szentgyorgyi & Furész 2007; Furész 2008) at- 
tached to the 1.5 m telescope at the Fred Lawrence 
Whipple Observatory on Mount Hopkins. The wave- 
length range 3900-9100 Å is covered in 51 orders at a 
resolving power of 44,000. The observations were made 
on November 2, 5, 11, and 12, 2020. Each spectrum 
is a combination of three 15-min. exposures. The flux 
for each exposure is about 400 photons per pixel. Small 
size of the input apertures (fibers) in TRES and CHI- 
RON rules out potential contamination from unresolved 
sources within a single TESS pixel. 


5.3. Spectral Analysis 


We have used the 12 spectra taken with CHIRON and 
TRES to extract the RVs of the primary stars in the 
three EBs of the system. The RVs of the sharp-lined 
primary of pair B are derived by the standard method, 
i.e., by cross-correlation of the spectra with a binary 
mask (CHIRON) or with a template (TRES) and fit- 
ting the resulting CCF. However, owing to the blending 
and rapid axial rotation, the CCFs are not suitable for 
measuring the RVs of the primaries in binaries A and C; 
instead, a different approach is needed based on model- 
ing the observed spectra. The light-curve analysis indi- 
cates that the secondary components in all three eclips- 
ing pairs are much fainter than the primaries. There- 
fore, we assume that the contribution of all secondaries 
to the spectrum is negligible and model it as a sum of 
three spectra of the primaries. The light-curve analysis 
indicates that the fluxes of all primaries in the TESS 
band are comparable. 

The RVs of the narrow B1 dip are determined by the 
standard procedure (approximation by a Gaussian func- 
tion). These RVs (as well as the 4 RVs from TRES) 
correspond to the circular orbit of B presented below. 

For modeling the spectra, we use the stellar parame- 
ters determined from the MCMC analysis (see Section 
6). Assuming synchronous stellar rotation with their re- 
spective orbits and zero obliquity, we compute equato- 
rial velocities of 48, 10.2, and 58 km s^!, which approx- 
imate the projected velocities because these binaries are 
eclipsing. 

The orders of the echelle spectra were merged together 
and normalized by the continuum. The merging proce- 
dure is not perfect and leaves residual waves in the con- 
tinuum in the area where the orders overlap. This minor 
defect is neglected here. The merged spectrum is con- 
structed on a logarithmic wavelength grid with a step of 
1 km s^! ranging from 480 nm to 650nm. The normal- 
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Figure 11. The observed spectra with a 10-pixel smoothing (red) are compared to their models (green) in the region around 
the Mg Ib triplet. Left: CHIRON spectrum, right: TRES spectrum. 


ized spectrum was also correlated with the same solar 
mask in the wider +400 km s^! range. These “wide” 
CCFs are slightly sub-optimal in comparison with the 
standard order-by-order CCFs in terms of the photon 
noise. However, the order-merged and normalized spec- 
trum is needed for modeling. The TRES spectra were 
transformed to the same order-merged format using the 
wavelength calibration provided in the headers. 

We use as a template the synthetic spectrum from 
the POLLUX library (Palacios et al. 2012) with Teg = 
6200K, logg = 4.0, and [Fe/H] = —0.5. The solar- 
metallicity template, chosen initially, has lines deeper 
than observed. The template is rotationally broadened 
using the calculated equatorial velocities. The instru- 
mental broadening is also included, but in the context 
of the present study it is negligible. Rotational broad- 
ening assumes a linear limb-darkening coefficient of 0.68 
(solar value). 

The orbital parameters of the primary stars in A and 
C were initially determined using the masses of the com- 
ponents from the MCMC system analysis (Section 6) 
and then iteratively improved; the orbit of B is well de- 
fined, and its RVs are assumed to be accurately known. 
Initially, we also assumed that the relative contributions 
of all stars to the spectrum are equal, but then refined 
the relative fluxes to A:B:C=0.3:0.4:0.3, making B the 
slightly brighter star. As we see below (Section 5.4), the 
two components resolved by speckle interferometry con- 
tain binaries AC and B, respectively. The magnitude 
difference in the J band of 0.27 mag, implies that the 
flux ratio AC:B=0.56:0.44, so B could be even a little 
brighter than assumed here. 

The fitting program selects one of the observed spec- 
tra and compares it to the model. The templates of the 
three stars are shifted by their respective RVs (known 
from the orbital elements) and by the barycentric cor- 


rection and summed in proportion to the assumed rel- 
ative fluxes. Figure 11 shows two examples comparing 
models to the observed spectrum. They match qualita- 
tively. Despite the smoothing, the observed spectra are 
noisy. A better assessment of the model is obtained by 
correlating it with the same solar mask and comparing 
the observed and modeled CCFs. This is illustrated in 
Figure 12 for two dates. The first spectrum, taken on 
JD 2459160, did not have a well-defined broad dip in 
the CCF because A1 and C1 had different RVs. On JD 
2459164, the dips of A1 and C1 overlapped, producing a 
clear signature in the CCF. Note that the dip amplitude 
of B1 apparently differs from the model by a variable 
factor. This could be caused by the fact that the star is 
a 0'4 visual binary, so the components can be mixed in 
slightly different proportions, depending on the guiding. 
'The CCF outside the dips is not constant; it varies ow- 
ing to random coincidences between spectral lines and 
mask. To some extent, this variation is captured by the 
modeled CCF. 

So far, the model uses the pre-computed RVs, without 
any fitting. Taking these RVs as the initial guess, we fit 
the RVs of A1 and C1 to minimize the sum of squares 
between the observed spectrum and its model. Fitting of 
the two parameters is done using the amoeba minimizer 
(Press et al. 1986). The RVs of B1 and other parameters 
(flux ratios, rotation speeds) are assumed known. A 
version of the code fitting all 3 RVs gives for B the same 
results as fitting the CCF dips. We also tried to fit 4 
parameters, including the relative flux of B. Its best-fit 
value of 0.36 is found consistently (the rms scatter is 
0.01) on all dates except 59164, when amoeba converged 
slowly and the best ratio returned by the code is 0.30. 

The RVs of A and C are used to determine their orbital 
parameters which, in turn, are used as the initial guess 
in further work. Depending on the details (initial guess, 
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Figure 12. Comparison between the CCFs computed for the real (green) and modeled (black) CHIRON spectra. 
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are computed independently, not fitted. The vertical lines mark the RVs of the components assumed in the model. 


fitting tolerance), the resulting RVs of A and C may 
differ by ~1 km s~!, except the dates where the dips 
of A and C strongly overlap and the differences may be 
larger. 

Minimization of the quadratic distance between the 
spectrum and its model is mathematically equivalent 
to maximizing their product, i.e., the cross-correlation. 
Owing to the artefacts of the spectrum and the intrin- 
sic mismatch between real spectra and templates, the 
residuals are much larger than the statistical errors. For 
the same reason, estimation of the errors of the derived 
RVs appears problematic. 

Table 3 lists the RVs of all 3 components derived from 
the CHIRON and TRES spectra. The RVs of B1 (mid- 
dle column) come from direct fitting of the CCF dip 
(CHIRON) or the CCF with a non-rotating template 
(TRES). The RVs of Al and C1 are determined by the 
spectrum modeling described above. The RVs of Al 
and C1 on BJD 59164, when their dips overlap, are less 
certain. 

The elements of circular spectroscopic orbits fitted to 
the RVs are also listed in Table 3, and the RV plots are 
shown in Figure 13. Each orbit is based on 12 RVs, 8 
from CHIRON and 4 from TRES. The TRES RVs are 
given a lower weight in the fits (with the latter rela- 
tive error bars taken to be 1.5 times larger than for the 
CHIRON points). The epoch To corresponds to the pri- 
mary eclipse, so the argument of periastron is fixed to 
w = 90?. An attempt to fit an eccentric orbit of B gives 
e = 0.005 + 0.005, so we assumed the orbit to be circular 
in subsequent analysis. 

Given the estimated masses of the primary-star com- 
ponents (see Section 6), the spectroscopic orbits con- 
strain the mass ratios. The final estimates of the com- 
ponents' masses are given in Section 6 using all available 


0.96 F =| 
0.94 F 7 
o92L JD 59164.838 I 
0.90 L. a cca [i aide ae ica ee ee es [ie Sg s 
—200 —100 O 100 200 
RV [km/s] 
Both CCFs 
Table 3. Radial Velocities 
BJD Al Bl C1 
— 2400000 (km s^!) 
CHIRON 
59160.7361 104.9 27.51 28.3 
59162.7152 84.1 28.30 107.5 
59164.7560 -16.8: 90.58 7.4: 
59165.7460 103.9 104.01 109.1 
59166.7443 43.0 94.01 140.5 
59167.7603 8.0 64.30 42.3 
59168.6616 114.9 34.78 1.9 
59175.7079 -2.8 73.15 87.9 
TRES 

59155.8896 67.7 70.84 6.7 
59158.8696 36.2 85.52 114.7 
59164.8419 -15.0 92.62 -7.6 
59165.8692 83.6 103.74 74.6 
Orbital Element 
Kı [km s^ !] 63.7 +1.5 44.28 +0.14 78.9 2.6 
Vo [km s7 +] 51.5 1.1 60.03+0.10 — 66.7 1.6 
To? 9151.868(6) 9151.446(4) 9150.193(7) 
a’ [km s !] 4 0.34 6 


Notes. (a) To is the epoch of the descending node of the RV 
curve (i.e., the eclipse time) in BJD - 2450000. (b) o is the 
rms residuals of the RV points from the fit. The rough 
estimated relative error bars on the individual RV points 
were scaled until X? per degree of freedom was unity. 


information, including the system SED and the analyses 
of the photometric lightcurves. 

Interestingly, the systemic velocities of A and C de- 
viate from the velocity of B in the opposite sense, and 
their mean, 59 km s^, is close to the velocity of B. This 
tells us that binaries A and C orbit each other with a pe- 
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Figure 13. Spectroscopic orbits. CHIRON RVs are plotted as black circles while the TRES RVs are indicated with green 


squares. 


Table 4. Measurements of AC,B at SOAR 


Date P.A. Sep. Am Filt. 

(JY) (deg)  (arcsec) (mag) 
2020.8236 257.74 0.4230 0.27 I 
2020.8368 257.61 0.4233 0.29 I 
2020.9243 257.62 0.4235 0.28 I 
2020.9243 257.69 0.4243 0.31 V 


riod of the order of several years, while binary B belongs 
to the visual secondary component (see Section 5.4). 


5.4. Speckle Imaging 


TIC 168789840 was observed with the speckle cam- 
era at the Southern Astrophysical Research Tele- 
scope (SOAR) on October 27, 2020 (JY 2020.8236). 
The instrument and data processing are described by 
Tokovinin (2018b). Several series of 400 images with 
exposure time of 24.4ms per frame were taken in the 
I band (824/170 nm) using the iXon-888 electron- 
multiplication CCD camera. The image cubes are pro- 
cessed by the standard speckle method. Orientation on 
the sky and pixel scale (15.81 mas) are determined from 
calibration binaries with well-known positions. The lat- 
est results from this instrument and references to other 
publications can be found in Tokovinin et al. (2020). 

The object was clearly resolved into a 042 pair at 
position angle of 257°7 with a magnitude difference 
AI = 0.27 mag (Figure 14). The true quadrant was de- 
termined from the shift-and-add images. Owing to the 
excellent 0"6 seeing on that night, the pair is partially re- 
solved even in the classical sense in the re-centered and 
co-added images produced from the data cubes. Ap- 
proximation of the semi-resolved classical image by two 
Moffat functions provides independent confirmation of 
the magnitude difference derived from speckle process- 
ing. Observation was repeated on November 1 and De- 
cember 3, 2020, and practically the same results were ob- 
tained (Table 4). Data over a wider field were also taken 
to ascertain the absence of other faint sources at larger 


separations, up to 8". The contrast limit for detection 
of other companions is about 4.0 mag at 1" separation 
and 5.5 mag at 3"and further out. 

One of the resolved components (the primary) is a 
close pair consisting of binaries A and C. However, 
separation of this inner pair should be less than ~30 
mas, otherwise it would be detectable by the asymme- 
try in the speckle power spectrum; no such asymmetry 
is found. 


6. ANALYSIS OF SYSTEM PARAMETERS 


The results of follow-up observations, along with the 
original T'ESS lightcurve and archival data, allowed for 
an extensive analysis of the system parameters. First, 
we determined a number of dimensionless ratios for each 
of the binaries, e.g., R/a and Teg ratios using PHOEBE 
(Prsa et al. 2011) and Lightcurvefactory (Borkovits 
et al. 2013; Rappaport et al. 2017; Borkovits et al. 
2018) from analysis of the disentangled photometric 
lightcurves. In a second step, we combined these ratios 
with the measured system spectral energy distribution 
(‘SED’) to determine the stellar parameters for all six 
stars using an MCMC analysis. These analyses are de- 
scribed in detail in the following sections. 


6.1. PHOEBE and Lightcurvefactory Analysis 


We start the analysis of the system parameters by first 
fitting the disentangled photometric lightcurves with 
two different binary lightcurve emulators: PHOEBE (Prsa 
et al. 2011) and Lightcurvefactory (Borkovits et al. 
2013; Rappaport et al. 2017; Borkovits et al. 2018). To 
further produce two independent sets of results, we use 
PHOEBE with the Fourier disentangled lightcurves (see 
Figure 5), and Lightcurvefactory with the iteratively 
disentangled lightcurves (see discussion in Section 3.1). 
In both of these analyses we used two simplifying as- 
sumptions: (i) circular orbits for all three pairs, and (ii) 
fixed effective temperatures for the primary stars in all 
three binaries. A logarithmic limb-darkening law was 
also applied for both analyses. 

The results of the PHOEBE fit to the Fourier dis- 
entangled lightcurves are shown in Figure 15, while 


Figure 14.  Speckle auto-correlation function of TIC 
168789840 (in negative rendering) recorded on 2020 Octo- 
ber 27 at SOAR. Two peaks B and B' on both sides of the 
center O indicate that it is a resolved pair; the true peak 
corresponding to the secondary component is marked by the 
white dot. The field size is 3/16, binary separation 07423; 
data taken in a wider field (up to 162) show absence of other 
fainter sources. The insert shows a long-exposure image pro- 
duced from the same data cube where the pair is partially 
resolved owing to the good 0"6 seeing. 


the fits to the iteratively disentangled lightcurves using 
Lightcurvefactory are shown in Figure 16. 

The resulting dimensionless parameters derived from 
these two fits are given in Table 5. These fits allowed 
for the determination of six values of scaled stellar radii, 
R/a (where a is the semi-major axis of the orbit), three 
values of primary to secondary Teg ratios, as well as the 
three orbital inclination angles. Furthermore, in order 
to get temperature ratios of the primaries of the differ- 
ent binaries with respect to each other, we made ad- 
ditional runs with Lightcurvefactory, simultaneously 
fitting the blends of any two of the three pairs (i.e. the 
time-series TAB, TAC, TBC - see Table 1). We regard 
the consistency between the two independent analyses 
of the dimensionless system parameters seen in Table 5 
to be quite encouraging. 
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Figure 15. PHOEBE fits of the fold of binary A (top), B 

(middle), and C (bottom). 


In addition, we find the ‘third light’ parameters for 
each binary, i. e., the amount of the extra flux contribu- 
tion over the flux of the eclipsing binary being consid- 
ered. Both PHOEBE and Lightcurvefactory have built- 
in functionality to solve for the third light parameter 
in any given lightcurve. We note, however, that the 
third light values are not particularly accurate at this 
stage of the analysis. This is remedied by the fact that 
the analysis described here, as well as the more com- 
plete MCMC analysis described in Section 6.2, are done 
iteratively, and the results become progressively more 
accurate upon iteration. 


6.2. MCMC Analysis of the Stellar Parameters 


We now combine the results of the dimensionless sys- 
tem parameters with several other pieces of information 
and constraints to solve for all of the stellar parameters 
for the six stars. Our approach is to fit for the six stel- 
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Table 5. Fitted Parameters Based on the T'ESS Photometric Lightcurves 


Fitted Parameter 


RVs/Iterative Disentanglement 


Lightcurvefactory PHOEBE 


Fourier Disentanglement 


Rai/aa 0.215 + 0.002 0.217 + 0.002 
Ra2/aA 0.093 + 0.004 0.087 + 0.003 
Rpi/ap 0.077 + 0.003 0.066 + 0.002 
Rp2/ap 0.031 + 0.002* 0.056 + 0.002* 
Rci/ac 0.233 + 0.009 0.246 + 0.007 
Re2/ac 0.094 + 0.011 0.107 + 0.005 
Teg, A2 / Test, A1 0.590 + 0.014 0.554 + 0.005 
Tas) Tees 0.692 + 0.009 0.681 + 0.004 
Teg,c2/Teg,ci 0.590 + 0.018 0.626 + 0.007 
Tog, A1/Teg, Bi 1.034 + 0.038 

Tog, A1 /Teg,ci 1.024 + 0.047 

Tacu Tae 0.953 + 0.037 » 
Inclination A [deg] 89.6 + 0.5 89.6+ 0.4 
Inclination B [deg] 88.5 + 0.5 88.1 + 0.3 
Inclination C [deg] 75.9+1 74.7 + 0.5 


Third Lights after Iteration! 


Third light to A 
Third light to B 
Third light to C 


0.707 + 0.038 
0.604 + 0.047 
0.688 + 0.057 


Notes. (a) The dimensionless quantities in this table for the three EBs were derived from the TESS lightcurves that were 
disentangled using two independent methods (see text for details), and two different binary lightcurve emulators 
(Lightcurvefactory and Phoebe). (*) The disparity in the radius of B2 from the two different approaches is due to the 
differences in the eclipse width and depth from their respective disentanglement methods. (+) The third light results are 
arrived at after iterating the lightcurve analysis as described in Section 6.1 with the MCMC analysis of the system parameters 
as described in 6.2 


lar masses and a common age, while making the explicit 
assumption that all the stars in the sextuple are coeval 
and that there has been no mass transfer among the 
constituent stars. We also employ as constraints (i) the 
measured SED for the system, (ii) MIST stellar evolu- 
tion tracks (Dotter 2016; Choi et al. 2016; Paxton et al. 
2011, 2015, 2019), and (iii) Castelli & Kurucz (2003) 
model atmospheres. When this analysis was carried out, 
there was no Gaia distance information for this object 
in DR2 (Lindegren et al. 2018). Therefore, we also fit 
for the distance to the source as well as the unknown 
interstellar extinction. 

In Table 6 we summarize exactly what the MCMC 
fitted parameters, the constraints, and the output pa- 
rameters are. In all, we are fitting 9 free parameters. 
On the other side of the ledger there are 12 easily iden- 
tified constraints (R/a and Tyg ratios, and RVs; see 
middle column of Table 6). In addition there are 26 


SED points, MIST evolution tracks”, Castelli & Kurucz 
(2003) model atmospheres?, and the assumption of a 
coeval evolution of the system without mass transfer. 
These latter items are hard to quantify in terms of a 
‘number’ of constraints; whether they are adequate will 
be determined by the uncertainties in the results. In the 
end, we hope to determine 21 independent parameters 
of the system, as listed in the 3rd column of Table 6. 
To carry out this fit for the 9 free parameters we used 
a Markov Chain Monte Carlo (MCMC, see, e.g., Ford 
2005) code modeled after the one used in Kurtz et al. 
(2020) and Rappaport, Kurtz, Handler, et al. (2020; 
submitted to MNRAS), but modified to handle six stars. 


2 The MIST tracks were for an assumed solar chemical composi- 


tion. 


3 The Castelli & Kurucz (2003) model atmospheres were also for an 
assumed solar chemical composition and for a fixed log g = 4.0, 
which well matches the primary stars in the problem (see Table 
7) that contribute 97% of the system light. 
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Figure 16. Lightcurvefactory fits of the fold of binary A 
(top), B (middle), and C (bottom). 


For the initial MCMC runs, the priors on the six stel- 
lar masses, the age and distance of the system, and the 
interstellar extinction were taken to be uniform over suf- 
ficiently large ranges so as to include all plausible values. 
For the final runs, the ranges of the priors were some- 
what narrowed, but the priors remained uniform over 
their respective ranges. In all cases, and for all param- 
eters, the range of priors was wider than +4ø of the 
finally determined parameter error bars. 

For each link in the MCMC chain we know the trial 
masses and the system age. From the evolution tracks 
we then also know all the corresponding radii and ef- 
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fective temperatures. The masses, combined with the 
known orbital periods, yield the semi-major axis of each 
of the three binaries. With this information we can 
check how well the R/a values and temperature ratios 
match the input values (see Section 6.1). The stellar 
radii and effective temperatures are then used in con- 
junction with the trial distance and Ay value, along with 
the atmosphere models, to compute the model compos- 
ite SED. These all contribute to x? in assessing whether 
to retry the step or make another jump from that point. 
For each addition to xy? we assume Gaussian distributed 
uncertainties in the RV values, the R/a values, temper- 
ature ratios, and the uncertainties in the SED points. 

We ran a dozen independent MCMC chains of 20 mil- 
lion links each to arrive at our results. Table 7 lists our 
fits to the system parameters, with uncertainties for the 
masses, radii, and Tugs. We also list in the Table sev- 
eral other parameters for each star that may be helpful 
in making sense of future RV or imaging observations 
of this system, e.g., the expected orbital velocities. Fig- 
ure 17 shows the posterior distributions for the six stellar 
masses, the six radii, and the six Tog values. The fourth 
panel in that figure gives the distributions of distance to 
the sextuple as well as of its age. 

The three short-period binaries would seem to be very 
similar ‘triplets’, each with a more massive primary (of 
~1.2 Mo) that is slightly evolved off the MS, and a 
secondary that is sub-solar and unevolved. The main 
difference among these three binaries is that one of them 
has an orbital period which is ~5 times longer than the 
other two. 

In Figure 18 we show the best fit to the SED data 
(from VizieR; Ochsenbein et al. 2000). The six thin 
curves are the contributions to the SED from the in- 
dividual stars. The heavy red curve is the sum of the 
contributions. The black points with error bars are the 
measured points. Note that the Galex (Morrissey et al. 
2007) NUV point is right on the model curve (though 
hard to notice). The best fit is for a distance of 571 pc 
in this figure (584 + 70 pc in Table 7) and an Ay value 
of 0.28. Now that the Gaia EDR3 (Lindegren et al. 
2020) are available, we have checked our fitted photo- 
metric distance with the parallax-determined value of 
593 + 150 pc, and find strong agreement. 

In Figure 19 we show the location of the six stars of 
'TIC 168789840 in the plane of stellar radius and effective 
temperature, with superposed stellar evolution tracks. 
All three of the primary stars lie close to the evolution 
track for a 1.2 Me star and have distinctly evolved away 
from the main sequence (between the TAMS and sub- 
giant phase). The three secondary stars are clearly sub- 
solar and near the main sequence. 
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Figure 17. MCMC outputs showing the distributions of system parameters. Note the similarity in each of the three primaries 
and secondaries, leading to our discussion of the binaries as “triplets”. 


To our knowledge, this is the first time that a fit for the 
properties of six stars, based largely on a set of overlap- 
ping photometric lightcurves and composite SED infor- 
mation, has been attempted. As a final demonstration of 
how the Lightcurvefactory photodynamical model us- 
ing all these parameters fits the original T'ESS lightcurve 
we show in Figure 20 the two curves superposed over a 
7-day segment of the TESS lightcurve. The correspon- 
dence with the actual data is quite gratifying. 


6.3. Inferences on the quadruple and sextuple orbits 


The SOAR speckle auto-correlation function (Fig- 
ure 14) shows two images of comparable brightness sep- 
arated on the sky by 07423, in agreement with, but more 
accurate than, the new Gaia EDR3 results (Lindegren 
et al. 2020) of 07374 + 07021. We had tentatively argued 
in Sections 5.3 and 5.4, that the brighter of the images 
was comprised of the A and C binaries, while the slightly 
fainter image (with AJ ~ 0.27 mag) was the B binary 
by itself. Here we further quantify that argument. 

We utilized the results of our MCMC analysis of the 
system parameters to predict the brightness of each bi- 
nary in the J band during each link in the MCMC chain. 
In Figure 21 we show distributions of the /-magnitude 
difference between the two images under the assumption 
that the inner quadruple is comprised of A+B, B+C, or 
A+C, respectively. In the first two cases, the measured 
SOAR magnitude difference of 0.27 mag has almost no 


plausible probability of agreement with model AT values 
of —0.94+0.31 mag and —0.88 4 0.17 mag, respectively 
(see Figure 21). By contrast, the hypothesis that A+C 
form the inner quadruple is consistent with the mea- 
sured value of AI with a value of —0.45 + 0.20 mag? 


The difference between the center-of-mass RVs of A 
and C also indicates that they are bound together in an 
orbit with a period of a few years. Thus, we are con- 
fident that the sextuple consists of an inner quadruple 
comprised of binaries A and C having a sky separation 
of S 30 mas, which in turn is orbited by binary B at a 
current-epoch sky projection of 07423. These two angu- 
lar separations amount to projected physical separations 
of S 18 AU and 250 AU, respectively. Circular orbits 
with these separations, coupled with the masses given 
in Table 7, would correspond to orbital periods of < 40 
yr and ~1700 yr, respectively. 

Depending on the exact separation of binaries A and 
C, there could well be observable Eclipse Timing Vari- 
ation (ETV) effects. We have therefore attempted an 
ETV analysis of the eclipse times of the binaries A and 


^ We have also verified that the G magnitude difference (Gac—Gp) 
reported in the new Gaia EDR3 release of -0.34, is in agreement 
with similar distributions from our MCMC analysis in the G 
band. 
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Figure 18. SED diagram for TIC 168789840. The six curves are the model contributions to the SED from the individual stars, 
while the heavy red curve is the sum of the contributions. The black points with error bars are the measured points (from 


VizieR; Ochsenbein et al. 2000). 
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Figure 19. The location of the six stars in TIC 168789840 in the plane of stellar radius and effective temperature. The tracks 
are taken from the MIST library (see Table 6 for references). The number next to each track is the corresponding stellar mass 


in units of Mo. 


C, similar to that used previously, e.g., in Zasche et al. periodic variation with a ~3.7 yr period, an amplitude 
(2019). of 0.0029 days, and an orbital eccentricity of 0.28. 

The top panel of Figure 22 shows the Observed minus The detection of a similar corresponding ETV for bi- 
Calculated (O — C) diagram for binary A, which has the nary C is tricky and yields rather uncertain results. The 
more readily detectable eclipses in the archival WASP reason is that the eclipses are too shallow and are barely 
and ASAS-SN data. Here there seems to be a fairly clear visible in the WASP and ASAS-SN lightcurves. The 


O — C diagram for binary C is shown in the bottom 
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Figure 20. Lightcurvefactory photodynamical model for the overlapping set of eclipses superposed on a 7-day segment of 
the T'ESS data. Blue points are the data in 30-min cadence while the red curve is the model. The different EB eclipses are 


marked (e.g. Ap marks the primary eclipse of binary A). 


panel of 22. Due to the relatively poor archival cov- 
erage of binary C, we adopted the following approach. 
Since we know the masses of both pairs A and C (Ta- 
ble 7), their respective amplitudes in the O — C diagram 
should follow the relation: a4/ac = Mc/Ma œ 10.14. 
Therefore, our joint analysis of both binaries used this 
simplification and both fits in Figure 22 were produced 
from a joint orbital solution, i.e., a fixed amplitude ratio, 
and common set of orbital parameters. As one can see, 
the predicted variation for pair C is difficult to discern. 
Much more precise times of eclipses, especially for pair 
C, are needed to confirm this hypothesis. Hopefully, new 
TESS data would serve as an ideal data source in this 
aspect. 

The orbital parameters we have found for the quadru- 
ple (AC) are given in Table 8. However, these should 
still be taken with caution especially due to poor archival 
data coverage of binary C. If we accept the ETV curve 
as the valid solution for the AC quadruple, then it di- 
rectly predicts the RVs of A and C as functions of time. 
In Figure 23 we show the expected RVs under the as- 
sumption that the systemic radial velocity, Vo, of binary 
B ~ 59 km s^! (see Table 3) also represents the center 
of mass velocity of the quadruple AC. We also plot on 
the figure the two values of systemic radial velocities, 
Vo, for binaries A and C (see Table 3) at the mean time 
the RVs were taken. One can readily see that there is a 
match for the predicted RVs if the orbital inclination of 
the AC quadruple orbit is of about 42°, in accord with 
what the ETV analysis indicates. 

For completeness, we note that it is highly unlikely 
that there are three unrelated EBs that are so precisely 


aligned along the line-of-sight just by chance. To calcu- 
late the probability corresponding to such a coincidence, 
we first compute the magnitudes of each EB from its flux 
contribution and find: 11.76, 11.58 and 11.98 Tmag for 
pairs A, B, and C, respectively. We then compare those 
with the number of nearby Gaia stars having magnitudes 
between 11.5 and 12.5 Tmag”, and with the speckle ob- 
servation from SOAR. There are 4 such stars in a 5’ x 5' 
region around the target: Gaia ID 4882948284462670000 
(Gmag = 12.11, Tmag = 11.68, using Stassun et al. 
2019; 4883001580713720000 (Gmag — 12.19, Tmag — 
11.76); 4882947498485530000 (Gmag — 12.77, Tmag — 
12.34); and 4883001615073460000 (Gmag — 12.8, Tmag 
= 12.37). Thus the probability of having one such 
star within 0.03" of the target (the SOAR limit on the 
separation of the inner quad)—and unrelated to it—is 
~ 4 x 107? and the probability for a star to be within 
4"4 of the target (the outer orbit separation as resolved 
by SOAR) is z 7 x 1076; thus the compound probability 
that TIC 168789840 is actually three unrelated EBs is 
z 3 x 10714. The equality between the RV of B and the 
mean RV of A and C is another strong argument that 
these stars are gravitationally bound in one system. 


7. SUMMARY 


In this work, we have presented the discovery of the 
first known sextuply-eclipsing sextuple star system TIC 
168789840. Our analysis shows that the orbital peri- 
ods of the three constituent eclipsing binaries are 1.570 


5 Assuming these are representative of the field of view. 


Table 6. System Parameters and Constraints in the MCMC 
Analysis 


Fitted Parameters Constraints“ Output 
Mai Rai/aa Mai 
Maz Ra2/aa Ma2 
Mpi Rgı/aB Mp1 
Mp2 Rg2/aB Mp2 
Moi Rci/ac Moi 
Mo» Re2/ac Mo» 


system age Test,A2/Teft,A1 system age 


distance Ter,p2/Tert,B1 distance 


extinction Ay Test,c2/Teft,c1 extinction Ay 


Kai Rai 
Kei Raz 
Ke Rei 
26 SED points Repo 
coeval assumption Re. 


MIST evolution tracks — Rc» 
Kurucz model spectra 


Te 
Te f,A2 
T. 


[o] 
mR 

wW 

= 


T. 
T. 


Teg,c2 


[o 
A 
w 
N 


[e 
R 

Q 

= 


Notes. (a) The R/a and temperature ratio constraints 
come from the light-curve emulator analysis of the 
disentangled TESS lightcurves for the three binaries. Kai 
is based on the RV analysis of the CHIRON and TRES 
spectra (see text). ‘Coeval’ assumption means that all six 
stars in the system are assumed to have been born at the 
same time, and that no mass transfer has occurred among 
them. The MIST stellar evolution models are from Dotter 
(2016); Choi et al. (2016); Paxton et al. (2011, 2015, 2019), 
while the ‘Kurucz’ model atmospheres are from Castelli & 
Kurucz (2003). 


days (binary A), 1.306 days (binary C) and 8.217 days 
(binary B), such that binaries A and C form an inner 
quadruple system with a period of about 4 years, and 
the latter forms the outer subsystem with a period of 
about 2,000 years. The three eclipsing binaries are prac- 
tically “triplets” with best-fit primary masses and radii 
of 1.23-1.30 Me and 1.46-1.69 Ro; secondary masses 
and radii of 0.56-0.66 Mo and 0.52-0.62 Ro; and pri- 
mary and secondary effective temperatures of 6350-6400 
K and 3923-4290 K, respectively. 

TIC 168789840 is a fascinating system that naturally 
merits additional observation and analysis. Though 
quite similar to the famous Castor system, the “triplet” 
nature of TIC 168789840 combined with the presence of 
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Figure 21. Brightness ratio distributions in each of the 
possible scenarios of the identity of the inner quadruple. For 
both the AB and BC quadruple possibilities, the measured 
value from speckle imaging lies on the extreme of the distri- 
butions. For AC, however, the measured value is well within 
the expected range. This provides strong supporting evi- 
dence to the RV analysis which concludes that AC is the 
inner quadruple. 


three primary and three secondary eclipses enable fur- 
ther investigations into its stellar formation and evolu- 
tion. Remarkable objects like TIC 168789840 or Castor 
give us insights on the formation of multiple systems — 
a matter of active research and debate. It is well known 
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Table 7. Computed Parameters for the Six Stars in TIC 168789840 


star Mass Radius Test Lumin a K vsini  logg 

(Mo) (Ro) K (Lo) (Ro) kms" kms" cgs 
A1 1.25 40.05 1.49+0.07 6400125 3.39 6.9 70.7 48.5 4.18 
A2 0.56 +0.04 0.52+0.04 3923+100 0.07 6.9 153.1 17.5 4.73 
B1 1.30 +0.08 1.69 34-0.22 6365+170 3.95 21.4 44.4 10.1 4.12 
B2 0.66 +0.03 0.62+0.02 4290+110 0.12 21.4 87.1 3.8 4.67 
C1 1.23 +0.10 1.45+0.28 6350+160 2.74 6.1 75.4 51.5 4.24 
C2 0.59+0.07 0.56+0.07 3990+190 0.07 6.1 154.3 20.9 4.72 
System dist age Av 

(pc) (Myr) (mag) 


584 + 70 3160 + 624 0.28 + 0.06 


a? 


Notes. All the parameters result from an MCMC study of the system constraints (see text). “a” is the binary semi-major axis. 


O-C diagram TIC 168789840 - pair A 


0.02 T T T T T T T T 
2006 2008 2010 2012 20147 2016 2018 2020 
0.01 + O T 
a D ü 
53 "MNI | NL 
| z * L i 
O +l NE | 
-001F L ee 
0.02 1 1 1 ae 1 1 1 
4000 5000 6000 7000 8000 9000 
BJD-2450000 
O-C diagram TIC 168789840 - pair C 
0.03 T T T T T T T T 
2006 2008 2010 2012 2014 2016 2018 2020 
0.02 - 
A 0.01 - 7 
a LT 
=: OF i dint 
O -0.01-f ^ 
-0.02 + 
0.03 


4000 5000 6000 7000 8000 9000 
BJD-2450000 


Figure 22. O — C diagram of binary A (top) and binary 
C (bottom). There are significantly more eclipse times for 
binary A than for C, and the former are more accurately 
measured. The orbits of the A and C binaries were fit jointly 
with a common set of orbital parameters (e.g., eccentricity, 
e, longitude of periastron, w, and time of periastron passage, 
7). The relative amplitudes of the two orbits were tied with 
a fixed ratio of a4/ac = Mc/MaA ~ 1. 


that components of hierarchical systems have correlated 
masses (Tokovinin 2018a), suggesting accretion from a 
common source. On the other hand, disk fragmenta- 
tion and subsequent migration, driven by accretion, ap- 
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Figure 23. Radial velocity predictions for the binary A and 
C center of masses based on the ETV solutions presented in 
Table 8 and Figure 22. The colored curves are for inclination 
angles of the AC quadruple of 10° (cyan), 30° (green), 50° 
(blue), 70° (red), and 90° (black), where the solid and dashed 
curves are for the C and A binaries, respectively. For an 
inclination of 42°, the expected values of Vo for the A and C 
binaries nicely match what one finds from the RV analyses 
of the A (solid circle) and C (open circle) binaries (see Table 
3). 


pears to be the dominant mechanism of close binary for- 
mation; its crude modeling can explain their statistics 
(Tokovinin & Moe 2020). The above model suggests a 
tight anti-correlation between the mass ratios of close bi- 
naries and the time of companion’s formation: low-mass 
companions are the latest to form. Formation of close 
binaries by several mechanisms acting separately or in 
combination and their migration can be ”observed” in 
numerical simulations of cluster collapse by Bate (2019). 

Transient massive disks prone to fragmentation likely 
result from an accretion burst, caused, e.g., by a close 
approach of two protostars surrounded by gas envelopes. 
Then one or both stars form secondaries via disk frag- 
mentation and also become bound together in a wide 
orbit (see Bate 2019). Accretion of the remaining gas 
drives the inner subsystems to closer orbits and shrinks 


Table 8. Fitted Parameters for the Inner Quadruple 


AC 
Parameter Value 
period 3.7 + 0.60 years 
Qa Sin? 0.516 + 0.110 au 
ac sint 0.510 + 0.110 au 
eccentricity 0.28 + 0.05 
WAC 166.0 + 25.2 
TAC 2457662 + 305 
f(M) 0.011 + 0.001 Mo 
i 42° 


Notes. The fitted projected semimajor axes, a, sin ? and 
ac sini were taken to have a fixed ratio in proportion to 
their measured inverse mass ratio. w and 7 are the 
argument of periastron and the time of periastron passage, 
respectively. f(M) is the mass function corresponding to 
the projected semimajor axes. The inclination angle, i, is 
inferred from the mass function and the measured masses 
of the A and C binaries. 


the outer orbit as well. Dissipative dynamics of an ac- 
creting triple or quadruple system evolving in a com- 
mon envelope presumably can align all its orbits in one 
plane. This scenario could explain formation of tight 
quadruples like VW LMi with nearly coplanar architec- 
ture (Pribulla et al. 2020) and compact coplanar triples. 

With regard to TIC 168789840, we might think that 
an encounter of the young binary AC with another star 
B led to its capture on a wide orbit, while strong accre- 
tion from the unified envelope, caused by this dynamical 
event, formed seed secondary companions to all stars by 
disk fragmentation. The seeds continued to grow and 
migrate inward, while the intermediate and outer orbits 
also evolved. The inner quadruple AC indeed resem- 
bles the tight coplanar quadruple VW LMi (outer pe- 
riod 1 yr), although in the latter the two inner mass 
ratios and periods (0.5 and 7.9 days) are not as simi- 
lar as in AC, and only one inner subsystem is eclipsing. 
This scenario, still speculative, explains the origin of the 
doubly-eclipsing inner quadruple AC and predicts that 
the orbit of AC should be coplanar with both inner bi- 
naries. The outer orbit of B around AC is much wider, 
and the eclipses of the binary B could be a matter of 
coincidence. In the Castor system, the outer orbit of 
^10 kyr period is not aligned with the 460-yr orbit of 
the intermediate quadruple, and only one of the three 
close binaries (the outer one) is eclipsing. 

For future measurements, we note that further reso- 
lution of the system may be possible with interferome- 
tetry. The axis of the inner quadruple AC is 4 AU or ~ 
7 mas, so might be marginally resolved by speckle at 10- 
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m telescope, or certainly at ELT. Consideration should 
also be given to GRAVITY, if fringe tracking is possi- 
ble. Future Gaia DR measurements have the potential 
to detect the ~ 4-yr wobble, with the DR3 catalog being 
released in December 2020. Additional RV monitoring 
will also give the AC spectroscopic orbit and upcoming 
TESS measurements may detect the ETV securely. 
Regarding our ongoing search for multiple star sys- 
tems, we continue to find more of these systems in the 
TESS data through a combination of machine learning 
(to limit the size of the data set) followed by a visual 
survey. T'ESS has allowed us to find well over 100 such 
candidate multi-star systems to date, with the analysis 
of another sextuple system by Jayaraman, Rappaport, 
Borkovits, Zasche et al. to follow this in this near future. 


NOTE ADDED IN MANUSCRIPT 


Since this paper was completed, we have received 
new TESS data which were taken at 2-minute cadence, 
thereby greatly improving the temporal resolution. We 
show in Figure 24 the new lightcurve and in Figure 
25 the folded, disentangled lightcurve. We have re- 
done those parts of the analyses which utilize the T'ESS 
lightcurve and find that the basic answers presented 
herein do not change significantly. 
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Figure 24. TESS lightcurve of TIC 168789840 in sector 31. 
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